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Theoretical design of global optimization algorithms can profitably utilize recent statistical me- 
chanical treatments of potential energy surfaces (PES's). Here we analyze a particular method to 
explain its success in locating global minima on surfaces with a multiple-funnel structure, where 
trapping in local minima with different morphologies is expected. We find that a key factor in 
overcoming trapping is the transformation applied to the PES which broadens the thermodynamic 
transitions. The global minimum then has a significant probability of occupation at temperatures 
where the free energy barriers between funnels are surmountable. 

PACS numbers: 61.46. +w,02.60.Pn,36.40.Ei 



Global optimization is a subject of intense current in- 
terest, since better optimization techniques can produce 
cost reductions and greater efficiency. It is therefore im- 
portant to understand the key requirements for a success- 
ful global optimization method, rather than proceeding 
on purely empirical or intuitive grounds. 

One of the difficulties associated with global optimiza- 
tion is the exponential increase in the search space with 
system size. For example, the number of possible con- 
formations of a typical protein is so large that it would 
take longer than the age of the universe to find the na- 
tive state if the conformations were sampled randomly 
(Levinthal's 'paradox' jlj]). This problem can be more 
rigorously defined using computational complexity the- 
ory; finding the global minimum of a protein [|2j or a 
cluster H is NP-hard. 

In fact, it is easy to design surfaces that result in ef- 
ficient relaxation to the global minimum, despite very 
large configurational spaces . Such surfaces are char- 
acterized by a single deep funnel j?J leading to the global 
minimum, and this feature may underlie the ability of a 
protein to fold to the native state 

Hence, global optimization methods should find sur- 
faces with a single funnel relatively easy to tackle. The 
problem is significantly harder if there is more than one 
funnel, since there are then competing relaxation path- 
ways. The extreme case would be when a funnel that 
does not lead to the global minimum is thermodynam- 
ically favoured down to low temperatures. On cooling 
the system would probably be trapped in the secondary 
funnel. This explains why naive simulated annealing will 
often fail. A number of atomic clusters bound by the 
Lennard-Jones (LJ) potential exhibit surfaces with just 
such a topography. In this paper we examine in detail 
one such case, an LJ cluster with 38 atoms, to show how a 
recently-developed global optimization method is able to 
overcome trapping. Our results indicate some necessary 
conditions the algorithm must satisfy if it is to succeed 
when applied to PES's with multiple funnels. 



Most small LJ clusters have global minima based upon 
Mackay icosahedra p0| . However, for 38 atoms the global 
minimum is a face-centred cubic (fee) truncated octahe- 
dron (Figure ), and for 75-77 and 102-104 the global min- 
ima are based on Marks' decahedra [jll]. These minima 
were initially discovered by construction |I^,[l3| , and un- 
til our recent application of a 'basin-hopping' algorithm 
Jl4j only the LJ38 global minimum had been found by an 
unbiased global optimization method. 
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FIG. 1. Energy profile of a pathway between the two low- 
est energy minima of LJ38, namely the fee truncated octa- 
hedron (bottom left) and a structure based on the Mackay 
icosahedron with C$ v point group symmetry (top right). 
2 1,/6 (T is the equilibrium pair separation of the LJ poten- 
tial. The method by which this pathway was obtained is 
described in Ref. Jl5[ . 

The multiple-funnel character of the LJ38 PES is re- 
vealed in the energy profile of a pathway between the fee 
global minimum and the lowest energy icosahedral mini- 
mum (Figure ) ]T5| . Maxima on this pathway correspond 
to true transition states (stationary points with a single 
negative Hessian eigenvalue), and the segments linking 
maxima and minima are approximate steepest descent 
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paths. The two lowest energy minima are well separated 
in configuration space, and to cross the barrier between 
them the cluster must pass through high energy minima 
associated with the liquid-like state of the cluster. Tran- 
sitions between the fee and icosahedral regions of configu- 
ration space can only occur if the high energy 'liquid-like' 
minima have a significant probability of being occupied. 
At low temperatures the Boltzmann weights for these in- 
termediate states are small leading to a large free energy 
barrier between the two regions. 
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FIG. 2. Equilibrium thermodynamic properties of the un- 
transformed LJ38 PES. (a) The probability of the cluster be- 
ing in the fee, icosahedral and 'liquid-like' regions of bound 
configuration space, (b) The heat capacity, C v . These re- 
sults were obtained by summing the anharmonic partition 
functions for a sample of minima appropriately weighted to 
compensate for the incompleteness of the sample Jl6[ . The 
liquid-like region of configuration space is defined as those 
minima with E > —171.6 e. 

Some of the equilibrium thermodynamic properties of 
LJ38 are shown in Figure . There is only a small en- 
ergy difference between the global minimum and the low- 
est energy icosahedral minimum. However, the entropy 
associated with the icosahedral structures is larger and 
icosahedra become favoured at low temperature. The 
centre of this transition is at a temperature of 0.12 efc -1 
(e is the pair well depth of the LJ potential and k is the 
Boltzmann constant); it gives rise to the small peak in 



the heat capacity (Figure b). The finite-size analogue of 
the bulk melting transition occurs at 0.18 e/c -1 , produc- 
ing the main peak in the heat capacity. These transitions 
are also reflected in the occupation probabilities for the 
different 'phases' (Figure a). 

At low temperatures one observes the cluster trapped 
in either the fee or icosahedral funnels, because of the 
large free energy barrier between them. On cooling from 
the liquid-like state, there is a thermodynamic driving 
force for entering the icosahedral region of configuration 
space. This effect is exacerbated by the topography of 
the PES. The structure of simple atomic liquids has sig- 
nificant polytetrahedral character 0Jl^l , and relaxation 
into the icosahedral funnel is more likely because icosa- 
hedra have more polytetrahedral character than the fee 
structures. Hence, the icosahedral funnel is directly con- 
nected to the low energy liquid- like states, whereas the 
fee funnel is only connected to the high energy liquid- 
like states. Therefore, it is extremely probable that the 
system will enter the icosahedral funnel on cooling and 
there become trapped. 

To observe the truncated octahedron we must simu- 
late the cluster in the small temperature window where 
both the high energy liquid-like minima and the fee struc- 
tures have small, but significant, probabilities of being 
occupied. Indeed, we did observe the truncated octahe- 
dron in this temperature range, but only once in regular 
quenches from a 0.25 /xs simulation for Ar. The situation 
for the larger clusters with non-icosahedral global minima 
is even worse. The decahedral to icosahedral transition is 
sharper and lies even further below the melting transition 
than for LJ 38 (e.g. for LJ 75 it occurs at T — 0.09 efc -1 ). 

The global optimization method that we analyze here 
belongs to the family of 'hypersurface deformation' meth- 
ods [Ts| . In this approach the energy function is trans- 
formed, usually to a smoother surface with fewer minima. 
The lowest minimum of this new surface is then mapped 
back to the original surface, but there is no guarantee 
that the global minima on the two surfaces are related 
and often there are good reasons to think that they are 
not ]I2| ]. In contrast, the transformation that we apply is 
guaranteed to preserve the global minimum. The trans- 
formed energy E is defined by: 



E(X) = min{E(X)} 



(1) 



where X represents the vector of nuclear coordinates and 
min signifies that an energy minimization is performed 
starting from X. 

The above approach is very similar to Li and Scher- 
aga's Monte Carlo (MC) minimization fl2(J. In our recent 
application to LJ clusters it has outperformed all other 
methods in the literature, finding all the known lowest 
energy LJ clusters up to 110 atoms, including those with 
non-icosahedral global minima Q. The method we use 
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to explore the E surface is simply a canonical MC simula- 
tion at constant temperature, in this case 0.8 ek" 1 . Inter- 
estingly, the other methods that have been most success- 
ful for LJ clusters are genetic algorithms |2l],^2| which use 
minimization to refine the new coordinates generated at 
each step, thus in effect performing a search on the same 
transformed surface, E. 
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FIG. 3. Equilibrium thermodynamic properties of the 
transformed LJ38 PES. (a) The probability of the cluster be- 
ing in the fee, icosahedral and 'liquid-like' regions of bound 
configuration space, (b) The configurational component of 
the heat capacity. 

The topography of the transformed surface is that of a 
multi-dimensional staircase with each step corresponding 
to the basin of attraction surrounding a minimum. The 
transformation has a significant effect on the dynamics. 
Not only are transitions to a lower energy minimum bar- 
rierless, but they can also occur at any point along the 
boundary between basins of attraction, whereas on the 
untransformed surface transitions can occur only when 
the system passes along the transition state valley. As 
a result intrawell vibrational motion is removed and the 
system can hop directly between minima at each step. 

However, the increase in interbasin transitions will 
not necessarily aid the location of the global minimum 
on a multiple-funnel PES unless the transformation also 
changes the thermodynamics. For LJ38 the probability of 



the system occupying the intermediate states between the 
funnels must be non-negligible under conditions where 
the icosahedral and fee structures also have significant 
occupation probabilities. This is indeed what is observed 
for the E surface (Figure ). The transitions have been 
smeared out and there is now only one broad peak in the 
heat capacity. Most significantly, the probability that 
the system is in the basin of attraction of the global 
minimum decays much more slowly, and the temperature 
range for which both the high energy liquid-like minima 
and the global minimum have significant probabilities 
is large. MC simulations anywhere in this temperature 
range easily locate the global minimum from a random 
starting point. Furthermore, the E transformation opens 
up new paths between the two lowest energy minima be- 
cause atoms are now able to pass through oneanothcr. 
Such paths would obviously not be feasible on the un- 
transformed PES. Consequently, even at T — 0.2 efc -1 
the cluster can escape from the icosahedral region of the 
transformed surface. 

We can understand the different thermodynamics for 
the two surfaces by writing the partition function as a 
sum over all the geometrically distinct minima on the 
PES. On the untransformed surface this gives for p s , the 
probability that the cluster is in minimum s, 



_ n s exp(-f3E s ) . ^ n s exp(-/3E s 

Ps\P) — -3JV-6 / 



v 



3N-6 



(2) 



where (3 = 1/fcT, N is the number of atoms, V s is the ge- 
ometric mean vibrational frequency of minimum s (cal- 
culated by diagonalizing the Hessian matrix), n s is the 
number of permutational isomers of minimum s and is 
given by n s — 2N\/h s , where h s is the order of the point 
group of minimum s. This equation is only approximate 
since it assumes that each well is harmonic. However, 
it does give a qualitatively correct picture of the ther- 
modynamics p3| , and although anharmonic corrections 
can be found ]16| they are rather cumbersome. For the 
transformed surface 

n s A s exp(-/3i? s ) 



Y, s n s A s exp(-f3E s 



(3) 



where A s is the hyperarea on the PES for which min- 
imization leads to minimum s. The A s values depend 
upon the available nuclear configuration space, which 
must be bounded to prevent evaporation. We have con- 
sidered both placing the cluster in a container and re- 
stricting the configuration space to hyperspheres around 
each local minimum by setting the coordinates to those 
of the relevant minimum in the Markov chain. Similar 
results are obtained for appropriate choices of container 
and hypersphere radii; the A s values reported below and 
the results in reference ]l4| ] were obtained by resetting 
the coordinates. 

Expressions (2) and (3) differ only in the vibrational 
frequency and A s terms. The fee to icosahedral and the 
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icosahedral to liquid transitions are caused by the greater 
number of minima associated with the higher tempera- 
ture state, which leads to a larger entropy. On the un- 
transformed surface this effect is reinforced by the de- 
crease in the mean vibrational frequencies with increasing 
potential energy, V fcc : v lcos : v Uquid = 1 : 0.968 : 0.864. 
Although these differences may seem small, when raised 
to the power 37V — 6 they increase the entropy of the 
higher energy states significantly, sharpening the transi- 
tions and lowering the temperature at which they occur. 
In contrast, A s decreases rapidly with increasing poten- 
tial energy, A fcc : A lcos : A Uquid = 1 : 0.0488 : 0.00122. 
These values were obtained by inversion of the occupa- 
tion probabilities obtained in MC runs of up to 10 6 steps. 
The decrease in A s reduces the entropy of the higher en- 
ergy states, causing the transitions to be broadened and 
the temperature at which they occur to increase. 

We can now explain in detail why the basin-hopping 
method is successful. First, the staircase transformation 
removes the barriers between minima on the PES without 
changing the identity of the global minimum, accelerating 
the dynamics. Second, it changes the thermodynamics of 
the surface, broadening the transitions so that even for 
a multiple- funnel surface such as that of LJ38, the global 
minimum has a significant probability of occupation at 
temperatures where the free energy barrier for passage 
between the funnels is surmountable. 

It is this latter feature which is especially important in 
overcoming the difficulties associated with multiple fun- 
nels and represents a new criterion for designing success- 
ful global optimization methods. 

The broadened transition that results from the trans- 
formation also differs markedly from the optimum condi- 
tions for protein folding, if we assume that proteins have 
evolved single-funnel surfaces in order to fold efficiently. 
A steep funnel provides a large thermodynamic driving 
force for relaxation to the global minimum , and also 
leads to a sharp thermodynamic transition. However, in 
global optimizations one cannot make assumptions about 
the topography of the PES and on a multiple-funnel sur- 
face such features can exacerbate the difficulties associ- 
ated with trapping in secondary funnels [||. 

Moreover, in protein folding there is the extra require- 
ment that the folded protein must remain localized in 
the native state, and this necessitates a sharp transition. 
There is no need for a global optimization method to 
mimic this property. Indeed, when applied to a PES 
with multiple funnels the optimum temperature for the 
basin-hopping approach is above the transition, where 
the system only spends a minority of its time in the global 
minimum. 

We have tested the basin-hopping method on a num- 
ber of other systems, and its performance is equally im- 
pressive. For example, it succeeds in finding the global 
minima for clusters bound by short range Morse poten- 
tials, which have much rougher energy landscapes than 



LJ clusters [Q. We have also obtained results for water, 
metal and silicon clusters which are helping to interpret 
experimental results. 
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